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Abstract The Aarhus code is the result of a long develop- 
ment, starting in 1974, and still ongoing. A novel feature is 
the integration of the computation of adiabatic oscillations 
for specified models as part of the code. It offers substantial 
flexibility in terms of microphysics and has been carefully 
tested for the computation of solar models. However, con- 
siderable development is still required in the treatment of 
nuclear reactions, diffusion and convective mixing. 
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1 Introduction 

What has become ASTEC started its development in Cam- 
bridge around 1974. The initial goal was to provide an im- 
proved equilibrium model for investigations of solar stabil- 
■ ity, fo llowing earlier work by IChristensen-Dalsgaard et al.] 
I ([l97J)- However, with the initial evidence for solar oscilla- 
tions and the prospects for helioseismology dChristensen-Dalsj 
l1976) the goals were soon extended to provide models for 
, comparison with the observed frequencies. Given the ex- 
pected accuracy of these frequencies, and the need to use 
them to uncover subtle features of the solar interior, more 
emphasis was placed on numerical accuracy than was per- 
haps common at the time. The code drew some inspiratio n 
from the Eggleton stellar evolution code (lEggletonLll97lh . 
which had been used previously, but the development was 
fully independent of that code. An early descrip tion of the 
code was given by IChristensen-DalsgaardI (Il982h. with fur- 
ther extensive details provided bv IChristensen-DalsgaardI 
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(Il978h : many aspects of this still stand and will only be sum- 
marized here. With the increasing quality and extent of the 
helioseismic data the code was further developed, to allow 
for more realistic physics; a major improvement was t he in- 
clusion of diffusion and settling (Christensen-Dalsgaar d et all 
il993h. This led to the so- called Model S of the Sun 
(IChristensen-Dalsgaard et afl, il996) which has found ex- 
tensive use as reference for helioseismic investigations and 
which at the time provided reasonably up-to-date represen- 
tations of the physics of the solar interior. In parallel, exten- 
sions have been made to the code to consider the evolution 
of stars other than the Sun; these include the treatment of 
convective cores and core overshoot, attempts to model red- 
giant evolution and the inclusion of helium burning. This 
development is still very much under way. 

For use in asteroseismic fitting a version of the code 
has also been developed in the form of a subroutine with a 
reasonably simple calling structure which also includes the 
computation of adiabatic oscillation frequencies as part of 
the computation. The combined package is available as a 
|ie t^i^ile^^iaking the installation relatively straightfor- 
ana ine"c6de has been successfully implemented on a 
variety of platforms. Nevertheless, it is sufficiently complex 
that a general release is probably not advisable. 



2 Equations and numerical scheme 

2.1 Formulation of the equations 

The equations of stellar structure and evolution are written 
with X = logjo q as independent variable, where logjQ is log- 
arithm to base 10; here q = in/M is the mass fraction where 
m is the mass interior to the point considered and M is the 
photospheric mass of the star, defining the photosphere as 
the point where the temperature is equal to the effective tem- 
perature. The equations are written on the form 
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Here r is distance to the centre, p is density, p is pressure, G 
is the gravitational constant, T is temperature, L is the lumi- 
nosity at r, £ is the energy generation rate per unit mass, H 
is enthalpy per unit mass, X]^ is the mass fraction of element 
k, Si]i is the rate of change of X]^ due to nuclear reactions, 
and '3lk and are the diffusion and settling coefficients for 
element ^. Also, the value of V = d log T/d log p depends on 
whether the region is convectively stable or not. In the case 
of a convectively stable region, 
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where a is the radiation density constant, c is the speed of 
light, and k is the opacity. The factor i// is discussed follow- 
ing Eq.[TT] below. The calculation of the temperature gradi- 
ent in convective regions is discussed in Section |4l In Eq.|5] 
the derivatives with respect to m should obviously be ex- 
pressed in terms of derivatives with respect to x. 

As discussed below, it is convenient to introduce the dif- 
fusion velocity on the form 



dm 



then Eqs|5]can be written as 
dXk ^Mq 
dx ~ &k 
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Here ^ = log 10, log being the natural logarithm. (In the 
code = M~^Yi( is used as variable.) For elements where 
diffusion and settling are ignored the original form, Eq. |5l 
with and "/^ set to zero is obviously used. 



2.2 Boundary conditions 

The centre, r = 0, is obviously a sing ular point of Eqs[T]-[4l 
As dis cussed in considerable detail by lChristensen-DalsgaardI 
this is treated by expanding the variables to second 
significant order in r to obtain the required conditions at the 
innermost point in the numerical solution. These are written 
in the form of expressions of the m and L at the innermost 
mesh point in terms of the other variables. The treatment 
of diffusion and settling has yet to be included consistently 
in this expansion; currently, conditions are imposed that es- 
sentially set the % to zero at the innermost meshpoint. The 



outer boundary of the solution of the full evolution equations 
is taken to be at the photosphere, defined as the point where 
T = Teff, the effective temperature. Consequently an obvious 
boundary condition is 



(10) 



where a is the Stefan-Boltzmann constant. To obtain a sec- 
ond condition, a expression is assumed for the dependence 
of T on optical depth T, conveniently written as 

r = 7;ff[T+^(T)]i/4. (11) 

Based on this the equation of hydrostatic equilibrium, on the 
form 



P = 
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dp^g_ 
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can be integrated assuming, currently, a constant gravita- 
tional acceleration g. The surface boundary condition 

2Tg 
K ' 

based on an approximate treatment of the outer parts of the 
atmosphere, is applied at the top of the atmosphere, T = Tmin, 
typically taken to be 10^^. The boundary condition is ob- 
tained by equating the pressure resulting from integrating 
Eq. [T2]to the pressure in the interior solution. 

The T{z) relation and Eq. [12] obviously define the tem- 
perature gradient V in the atmosphere . To ensure a contin - 
uous transition to the interior 1 follow iHenyey et all (Il965h 
and include the factor \j/ = 1 +^'(t) in Eq.|6l 

2.3 Numerical scheme 

As indicated in Eqs[T]-|5]and|8] |9] the dependent variables 
on the left-hand sides of the equations are expressed in terms 
of the set 



(9) {y,} = {logio'',logioP,,logio7',logioi,^/i,>i:'} , 



i = l,...,I. 



(14) 



Here the index k runs over all the abundances considered, 
whereas the index k' runs over those elements for which 
diffusion is included. In the solution of the equations, the 
computation of the right-hand sides is the heaviest part of 
the calculation and hence needs to be optimized in terms 
of efficiency. As discussed below, this may involve express- 
ing the thermodynamic state in terms of a different set of 
variables {z/}, related to the {y,} by a non-singular trans- 
formation. Using also the reformulation, Eqs [8] and |9] of 
the diffusion equation, the combined set of stellar-evolution 
equations consists of equations of one of the following three 
types: 
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u = h+h + \,...h+h + h 



Type III 



(17) scheme is used, e.g., for the time derivatives in the energy 
equation where the characteristic timescale is much shorter 

(18) than the evolution time and hence short compared with the 

(19) typical step At in time. 



The equations of Type 111 are obviously relevant to the evo- 
lution of abundances of elements for which diffusion is not 
taken into account. To this must be added the transformation 



yi=yi{x;zj;t) , «'=1, 



(20) 



In practice, z, = for / ^ 2; the choice of Z2 depends on 
the specific equation of state considered. The equations are 
solved on the interval [xi ,X2], with boundary conditions that 
can be expressed as 



ga{x\;zi{x\)) = , a = l,...,KA, 



(21) 
(22) 



The equations are solved by means of the Newton- 
Raph son-Kantorovich scheme (see lRichtmyeiifl957l : lHenriciI 
in the s tellar-evolution context kno wn as the Henyey 
scheme (e.g., iHenvev et al.L Il959l 1 19641) . We introduce a 
mesh jci = < x^ < . . .x^ = X2 in x and consider two 
timesteps f ' and f '"'"', where the solution is assumed to be 
known at timestep t\ Also, we introduce 
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with a similar notation for A, 
timestep r'"'"^ Then Eqs[l5]-[ 
ing difference equations: 
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as well as for quantities at 
I are replaced by the follow- 
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-yr = 4f^[0„/r+' + (l-0„)/r] , 



n = l,...,N; u = h+l2 + \,...h+h+h ■ (26) 

Here Ax" = x"+^ -x" and At' = t'+^ - t\ Also, the pa- 
rameters di allow setting the centralization of the difference 
scheme in time, 0, = 1/2 corresponding to time-centred dif- 
ferences and di = 1 to a fully implicit scheme. The former 
clearly provides higher precision but potentially less stabil- 
ity than the latter; thus time-centred differences are typically 
used for processes occurring on a slow timescale, such as 
the change in the hydrogen abundance, whereas the implicit 



together with the 
are solved using 



The algebraic equations |24] - [ 
boundary conditions, Eqs |2T] and 
Newton-Raphson iteration, t o determine the solution jz"'''^' } 
at the new time step (see also IChristensen-Dals gaardl 1 1 982l) . 
Given a trial solution {zj'^^^} the equations are linearized in 

terms of corrections 5z"'^~^^ = z" - z" and the result- 
ing linear equations are solve d using forwards elimination 
and backsubstitution (e.g.. Ba ker et al.L [1971). This process 
is repeated with the thus corrected solution as trial, until con- 
vergence. 

The initial ZAMS model is assumed to be static and with 
a prescribed chemical composition. Thus in this case only 
Eqs|24]are solved. Time evolution is started with a very small 
At for the initial non-zero timestep. 

The number of meshpoints is kept fixed during the 
evolution, but the distribution of points is varied accord- 
ing to the change in the structure of the model. The mesh 
algorithm is based on t he firs t-derivativ e stretchiiig intro - 
duced by iGough et all ( Il975h (see also lE ggletonl [l97lh . 
taking into account the variation of several relevant vari- 
ables. An early version of the implementation was described 
by Christensen- Dalsgaard ( 1982), but this has subsequently 
been substantially extended and is still under development. 
In particular, a dense distribution of points is set up near the 
boundaries of convection zones, although so far points are 
not adjusted so as to be exactly at the edges of the zones. 
After completing the solution at each timestep the new mesh 
is determined and the model transferred to this mesh using, 
in general, third-order interpolation; linear interpolation is 
used near boundaries of convective regions and in other re- 
gions where the variation of the variables is not smooth. 

Having computed model at timestep f^+^ the next time- 
step is determined from the change in the model between 
timesteps and f*+'; this involves a fairly complex algo- 
rithm limiting the change in, e.g., logjg/7, logjgT and X 
at fixed m. To compensate for the fairly crude treatment of 
the composition in a possible convective core, the change 
of X in such a core is given higher weight than the gen- 
eral change in X. The algorithm correctly ensures that very 
short steps in time are taken in rapid evolutionary phases. 
In typical simple calculations, assuming ^He to be in nu- 
clear equilibrium, roughly 35 (100) timesteps are needed to 
reach the end of central hydrogen burning in models with- 
out (with) a convective core, and 100 (200) steps to reach 
the base of the red-giant branch. Calculations with more 
complex physics or requiring higher numerical precision ob- 
viously require a substantially higher number of timesteps. 
Evolution up the red giant branch typically requires a large 
number of timesteps owing to the rapid changes at fixed m in 
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the hydrogen-burning shelQ ahhough the timestep algorithm 
has options to reduce the weight given to this region. 



3 Microphysics 

As the code has developed over the years a number of op- 
tions have been included for the microphysics, although not 
all of these have been kept up to date or properly verified. 
As a general principle, the code has been written in a mod- 
ular way, so that replacing, for example, routines for equa- 
tion of state or opacity has been relatively straightforward. It 
should be noted that the use of a different set {z,} of de- 
pendent variables on the right-hand side of the equations 
yields addit ional flexibility in the repl acement of aspects of 
the physics. iDappen and GuzikI (l2000l) provided a review of 
the treatment of the equation of state and opacity in stellar 
modelling. A review of nuclear reactions in the solar interior, 
of relevance in the rnore g eneral stellar case, was provided 
by lParker and Rolfsl (Il99lh . 



3.1 Equation of state 

The o riginal version of the code used the lEggleton et al.l 
(Il973h equation of state and that remains an often used op- 
tion. In this case zi = logio/^ where /, introduced by Eggle- 
ton et al. and related to the electron-degeneracy parameter, 
is used as one of the thermodynamic variables; this allows 
explicit calculation of partial ionization and hence a very 
efficient evaluation of the required thermodynamic quan- 
tities. The formulation also includes a crude but thermo- 
dynamically consistent implementation of 'pressure ioniza- 
tion' (which actually results mainly as a result of the high 
density in deep stellar interiors) which provides apparently 
reasonable results. The fact that the whole calculation is 
done explicitly makes it entirely feasible, if somewhat cum- 
bersome, to evaluate analytical derivatives. Up to second 
derivatives of pressure, density and enthalpy are provided 
in a fully consistent manner, whereas third derivatives, re- 
quired for the central boundary condition, assume full ion- 
ization. Partial electron degeneracy is included in the form 
of expansions that cover all levels of degeneracy and rela- 
tivistic effects. 



However, the resulting equation of state captures major as- 
pects of the more complex forms discussed below. 

More realistic descriptions of the equation of state re- 
quire computations that are currently too complex to be in- 
cluded directly in stellar evolution calculations. Thus inter- 
polation in pre-computed tables is required. The first such 
set to be included wa s the so-called MHD equation of state 
(iMihalas etalill990h . based on the chemical picture where 
the thermodynamic state is obtained through minimization 
of an expression for the free energy including a number 
of contributions. Application of this formulation to a com- 
parison with observed solar oscillation frequencies showed 
a very substantial imp rovement in the agreement betwee n 
the Sun and the model d Christensen-Dalsgaard et allll988[) . 
Further updates to the MHD equation of state have been 
made but they have so far not been implemented in ASTEC. 
An alternative des cription is provided b y the so-called OPAL 
equation of state dRogers et al.l 1 19961) . based on the physi- 
cal picture where the thermodynamic state is obtained from 
an activity expansion. This has been the preferred equa- 
tion of state for solar modelling with ASTEC, used, e.g., in 
Model S. 

The early versions of both the MHD and OPAL equa- 
tions of state suffered from a negl ect of relativistic effects 
in th e treatment of the electrons dEllio tt and KosovichevL 
19981). This has since been corrected (Gong et all bOOll : 
Rogers and Navfonovl |2002|) . Also, the OPAL formula- 
tion suffered from inco nsistencies between some of the 
variables prov i ded ( e.g. iBoothrovd and Sackmannl I2003t 
'Scufla ire et al.L l2007h : effects of the inconsistency on the 
computation of adiabatic oscillations are discussed bv lMova et al.l 
(f2007h . This has been improved in the latest version of the 
OPAL equation of state0 which has also been implemented 
in ASTEC. 

Interpolation in the OPAL tables uses quadratic interpo- 
lation in p, T and X. Typically, a single representative value 
of Z is used, even in cases with diffusion and settling of 
heavy elements, although the code has the option of using 
linear interpolation between two sets of tables with different 
Z. 



3.2 Opacity 



As an upgrade to the basic EFF formiil^tinn lchristpnspn-naH^aeiiigMyH t ij^sJma nt of th e opacities used tab les from 

ICbx and Stewarti (119701 and Cox and Tabor ( 1976) with an 
interpolation scheme developed by Cline (1974). A major 
improveme nt was the implementatio i i of the OPAL opac - 
ity tables ( Rogers and Iglesiasl 1 19921 : llglesias et all 1 19921) . 
With various updates of the tables this has since been the 
basis for the opacity calculation in the code. The most re- 
cent tables, including a variety of mixtu res of the heavy ele- 
ments are based on the computations bv lRogers and Iglesiasl 
(Il995h . Atmo spheric opacit i es must be supplied s eparately; 
here tables by lKuruczl d 1 99 ih . [Alexander and Ferguson (1994h 

^ see 

|http: //www-phys . llnl .gov/Research/OPAL/opal .htmll 



(fl992l) included Co ulomb effects, in the Debye-Huckel for 
mulation, following iMihalas et al.l (n"988i) ; unlike some ear- 
lier treatments of the Coulomb effect this ensured that ther- 
modynamic consistency was maintained; as a result, a sub- 
stantial effect was found also on the degrees of ionization 
of hydrogen and helium, resulting from the change in the 
chemical potential. The computation of the Coulomb effects 
consequently requires some iteration, even if the EFF vari- 
ables are used, and hence some increase in computing time. 

' This problem is avoided in the implementation by Eggleton ( 1 971h 
where the equations are solved using an independent variable that is 
directly tied to such strong variations in the model. 
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or 



Her guson etall (l2005h have been used, with a smooth and 



matching to the interior tables at logjg T = 4. Ele ctron con- 
ductio n may be included based on the tables of lltoh et"al] 
(Il983h . 

Interpolation in the OPAL tables i s carried out with 
schemes de veloped by G. Houdek (see iHoudek and Rogll 
Il99lll996h . The tables are provided as functions of {R, T), 
where R = p/T^. For interpolation in (log/?, log T) two 
schemes are available. One uses a minimum-norm algorithm 
with interpolating function defined in a pi ecewise f ashion 
over triangles in the log/? — log T plane dNielso n'. '1980, 



1983 ). The second scheme uses birational splines (Spath, 



19911) . In practice the latter scheme has been used for most of 



the model calculations with ASTEC. Interpolation i n X and 
logZ is carried out using the univariate scheme of I Akimal 
il97 (f). Note that the relative composition of Z is assumed 
to be fixed; thus differential settling or changes in heavy- 
element composition resulting from nuclear burning are not 
taken into account. 

The code includes flexible ways of modifying the opac- 
ity, to allow tests of the sensitivity of the model to such mod- 
ifications. An extensive survey of the response of solar mod- 
els to localized opacit y changes, specified as functions of 
login r , was made by iTripathy and Christensen-DalsgaardI 
(119981) . 



3.3 Nuclear reactions 



The calculation of the nuclear reaction rates is based on the 



are included; these play an important role in ensuring the 
gradual conversion of '^O to I'^N and hence increasing the 
importance of the CN cycle. This part of the cycle is as- 
sumed to be fully controlled by the '^0('H,7)'^F reaction 
and the branching ratio between the ^^N('H,'^He)'^C and 
'5N(iH,y) 1^0 reactions. 

Helium burn ing has been included in the code using the 
expressions of Angulo et al.l (IT999) . and including also the 
reaction '^C(^He, 7)^^0. However, the code is unable to deal 
with helium ignition in a helium flash. Thus models with 
helium burning are restricted to masses in excess of 2.3M0 
where ignition takes place in a relatively quiet manner. Also, 
as in cases of hydrogen burning (cf. Sect. |4|, the treatment 
of semiconvection that may be required in this phase has not 
been implemented. 



3.4 Diffusion and settling 

Diffusio n and settling are treated in th e approximations pro- 
posed bv lMichaud and ProffittI (Il993[) . with revisions kindly 
provided by Proffitt. For completeness I give the complete 
expressions in Appendix A. If included, diffusion of heavy 
elements assumes that all elements behave as fully ion- 
ized '^O; this is a reasonable approximation in the solar 
case where the outer convection zone is relatively deep, but 
becomes increasingly questionable in more massive main- 
sequence stars. Here, also, effects of sel ective radiative levi- 



usual approximations to the reaction integrals (e.g.,'Clavton', tation should b e taken into account (e.g. . [Richer et al.L[l998t 
1968) u sing a variety of c oeffic i ents f Bahcall and Pinsonneau ljjTurcotte et al.L [l998); there are no current plans to do so in 



19951 ; lAdelberger et all 1 19981 ; |Angulo et al., 1999). Elec- 



tron screening is computed in the Salpeter (1954) approx- 
imation. Electron ca pture by ^Be is treated according to 
iBahcall and Moelled dl969t) . 

The nuclear network is relatively limited and is one of 
the points where the code needs improvement. This is to 
some extent a heritage of its origin as a solar-modelling 
code, as well as a consequence of the fact that pre-main- 
sequence evolution is not computed. In the pp chains ^D, 
^Li and ^Be are assumed to be in nuclear equilibrium. On 
the other hand, the code has the option of following the 
evolution of ^He, although in many calculations it is suffi- 
cient to assume ^He to be in nuclear equilibrium. To simu- 
late the evolution of the ^He abundance X(^He) during the 
pre-main-sequence phase the initial zero-age main-sequence 
mode assumes the X(^He) that would have resulted from 
evolution over a specified period ?3h„ at constant co nditions, 
as described by lChristensen-Dalsgaard et"an d 19741) but gen- 
eralized to allow a non-zero initial abundance. A typical 
value is fsfjg = 5 X 10^ yr. 

In the CNO cycle the CN part is assumed to be in nuclear 
equilibrium and controlled by the rate of the ^'^N('H, y)'^0 
reaction. In addition, the reactions 
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0(^H,7)i^F(e+Ve)"0(^H,'^He)>'*N 



the code. Various approximati ons to turbulent diffusion ca n 
be included, partly inspired bv lProffitt and MichaudI (Il99lh . 
In addition, the code has the option to suppress settling in the 
outer parts of the star, to allow modelling of diffusion and 
settling in the cores of relatively massive stars where other- 
wise settling beneath the thin outer convection zone would 
result in a complete depletion of the surface layers of ele- 
ments heavier than hydrogen. 

At present, diffusion and settling is coupled to nuclear 
evolution in the consistent manner of Eq.|5]only for helium. 
For the remaining elements taking part in the nuclear net- 
work diffusion is neglected. Correcting this deficiency is an 
obvious priority. 



4 Treatment of convection 

The temperature gradient in con vect ion zones is usually 
computed using the Vitens^ dl953h and lBohm-Vit ense ( 1958j) 
version of mixing-length theory; for completeness, the ex- 
pressions used are provided in Appendix B. The mixing 
length is a constant multiple, ttuhUp, of the pressure scale 
height Hp. In addition, emulations of the Canuto and Mazzitellil 
( fl991,) formulation, established by .Monteiro et al. (1996) . 
can be used. 
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Convective regions can obviously, at least in stars that 
are not extremely evolved, be assumed to have uniform com- 
position. This can in principle be achieved by including a 
very high diffusion coefficient in such regions. In ASTEC 
this is used in convective envelopes, ensuring that they are 
chemically uniform. The treatment of convective cores re- 
mains a concern and an area of active development, how- 
ever. Given the lack of a proper treatment of the diffusion of 
all elements an explicit calculation of the chemical evolution 
is required. This is characterized by the (assumed homoge- 
neous) abundances Z<. ^ of the elements. The rate of change 
of these abundances can be written 

^ = J', -f — ^ [X,{x,,) - X,, e] + — ?^(xcc) ; (27) 

Clf ^cc ^cc 

here q^c is the mass fraction in the convective core, Xqc = 
logio(?cc), and 

^k = —r^^kdq (28) 
qcc Jq 

is the reaction rate averaged over the convective core. In the 
second term in Eq. |27]Xi(ji:cc) is evaluated just outside the 
core; this term only has an effect if there is a composition 
discontinuity at the edge of the core, i.e., if the core is grow- 
ing and diffusion is neglected. Finally, the term in fit(jCcc) 
is obviously only included in cases where diffusion is taken 
into account. 

In models with a growing convective core, and neglect- 
ing diffusion, a discontinuity is set up in the hydrogen abun- 
dance at the edge of the core (see also Fig. [Hi. This situ- 
ation arises in intermediate-mass stars (with masses up to 
around 1.7Mq) where the gradual conversion during evolu- 
tion of '^O into '^N causes an increase in the importance 
of the CNO cycle in the energy generation (for more mas- 
sive stars the CNO cycle dominates even with the initial '^N 
abundance). Since pressure and temperature are obviously 
continuous, there is also a discontinuity in density and opac- 
ity K, leading also to a jump in the radiative temperature 
gradient V^ad (cf. Eq. |6ll. Since k increases with X and p, 
while p decreases with increasing X, it is not a priori clear 
how K and hence Viad react at the discontinuity; in practice 
the common behaviour is that Vjaj increases going from the 
value of X in the convective core to the higher value in the 
radiative region just outside it. This raises the question of 
the definition of convective instability: if the border of the 
convective core is defined using the composition of the core, 
the region immediately outside the core is therefore convec- 
tively unstable. As a consequence, ASTEC defines the bor- 
der of the convective core by the abundance in the radiative 
region, leaving a small convectively stable region just be- 
low the border, which nevertheless is assumed to be fully 
mixed in the standard ASTEC implementation. This may 
be regarded as an example of semicon vection, of somew hat 
uncertain physical consequences (e.g.. iMerrvfieldl [19951) . A 
different scheme, now under implementation, is discussed in 
Section [6] 

Various options exist for the treatment of convective 
overshoot. Simple formulations involve compositionally 



fully mixed overshoot regions from the convective core or 
convective envelope, over a distance (XovHp below a convec- 
tive envelope, or aovniin(rcc,/^p) above a convective core, 
where Tec is the radius of the core. The overshoot region may 
be taken to be either adiabatically or radiatively stratified. 
A more elaborate treatment of overshoot fro m a convective 
envelo pe has been implemented, following iMonteiro et alj 
(119941) . where various forms of the temperature gradient can 
be specified, still assuming the overshoot region to be fully 
mixed. Thi s is being exten ded to emulate the overshoot sim- 
ulations bv lRempell (120041) . 



5 Implementation details 

When computing models of the present Sun it is important 
to adjust the input parameters such as to obtain a model that 
precisely matches the observed solar radius, luminosity and 
ratio Zs/Xs between the abundances of heavy elements and 
hydrogen, at the present age of the Sun. This is achieved 
by adjusting the initial hydrogen and heavy-element abun- 
dances Xq and Zo as well as the mixing-length parameter 
O'ml (or another parameter characterizing the treatment of 
convection). In ASTEC the iteration to determine these pa- 
rameters is carried out automatically, making the computa- 
tion of solar models, and the exploration of the consequences 
of modifications to the input physics, rather convenient. 

The ASTEC code has grown over three decades, with a 
substantial number of different uses along the way. This is 
clearly reflected in the structure of the code as well as in 
the large number of input parameters and flags that control 
its different options. These are provided in an input file, in 
many cases using simply the defaults provided in the source 
of the code. Also, several different executables can be pro- 
duced, reflecting partly the evolution of the code and partly 
different versions of the physics, in particular the equation 
of state and the opacity, as well as the option to include 
diffusion and settling. To make the code somewhat more 
user-friendly, scripts have been made which allow simply to 
change a few key parameters, such as mass, heavy-element 
abundance and number of timesteps, by editing templates of 
the input files. Consequently, the code has been used with 
success by several users, including students and postdocs at 
the University of Aarhus. 

In addition to summary output files listing global prop- 
erties of the models in the evolution sequence, output of de- 
tailed models, on the full mesh of the calculation, can be 
made in three different forms: the so-called emdl files, in- 
cluding just the variables {z,} actually used in the calcula- 
tion, as well as a complete listing of the input parameters, to 
provide full documentation of the calculation; the sundl files 
which provide the va riables needed for the Aarhus adiabatic 
pulsation code (see .Christensen-Dalsgaardl 12007 al) ; and the 
gong files, giving an extensive set of variables for use in fur- 
ther investigations of the models or plotting. A convenient 
way to use the code, without overloading storage systems 
with the large gong files, is to store the full emdl file and 
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subsequently read in models at selected timesteps, to output 
the corresponding amdl or gong files. 

For asteroseismology it is evidently crucial to compute 
oscillation frequencies of the computed models. The full cal- 
culation of frequencies, for given input parameters to the 
evolution calculation, often needs to be carried out as part 
of a larger computation, e.g., when fitting observed frequen- 
cies to determine the properties of the model. To facilitate 
this a version of the code has been made where the evolution 
calculation and all parts of the adiabatic oscillation calcula- 
tion is carried out by a single subroutine call, with internal 
passing of the intermediate products of the calculation. This 
subroutine can then be called by, for example, a fitting code. 
An example of such use is the applicati on of the code in 
genet ic-algorithm fitting (e.g., Metcalfe and CharbonneauL 
l2003h . under development by T. Metcalfe, High Altitude Ob- 
servatory. 

The code has been implemented on a variety of platforms 
and appears to be relatively robust. To simplify the instal- 
lation, a complete tar package including all the required 
files, with a setup script and a full makefile, has been estab- 
lished. However, the complexity of the code and the lack of 
adequate documentation makes it unrealistic to release it for 
general use. 



6 Further developments 

From the preceding presentation it is obvious that there are 
significant deficiencies in ASTEC, and work is ongoing to 
correct them. A fairly trivial issue is the restricted nuclear 
network and the failure to include all elements undergoing 
nuclear reactions in the full diffusive treatment. Rather more 
serious problems concern the treatment of the borders of 
convective regions; even though this obviously also involves 
open issues of a basic physical nature, the code should at 
least aim at treating these regions in a numerically consis- 
tent, even if perhaps not physically adequate, manner. A 
serious problem is the failure of the code for models with 
convective cores, when diffusion an d settling of both helium 
and he avy elements are included (cf. lChristensen-Dalsgaardl 
l2007bh : on the other hand, the case of just helium diffusion 
can be handled. This problem may be related to issues of 
semiconvection where convective stability is closely related 
to the details of the composition profile dMontalban et all 
l2007h . 

In ASTEC I have considered two cases of semiconvec- 
tion, although the implementation is still under develop- 
ment. One case concerns growing convective cores in non- 
diffusing models (see also Sect. |4]l, as illustrated in Fig. [T| 
The dashed curves illustrate the properties of a model com- 
puted with the normal ASTEC implementation, where the 
hydrogen abundance X is discontinuous at the boundary of 
the core (see panels a and b). As discussed above, the extent 
of the convective core is determined by the behaviour of the 
radiative temperature gradient V^ad in the radiative region 
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Fig. 1 Hydrogen-abundance profile and Vrnd — against fractional 
radius in two 1.5M .j models of age 1.36 Gyr and with Z = 0.02. 
The dashed curves illustrate the usual treatment in ASTEC, where the 
boundary of the convective core is determined by the composition in 
the radiative region just outside it, and the hydrogen abundance X is 
discontinuous. In the model illustrated by the solid curves a gradient in 
X is set up such that Vjad — ^ad = in the outer parts of the convective 
core. 

(see panel c)|3 leading to a convectively stable region just be- 
neath the boundary which nonetheless is assumed to be fully 
mixed. A possibly more reasonable treatment, illustrated by 
the continuous curves, assumes that a hydrogen-abundance 
profile is established such that Vjad = ^ad in the outermost 
parts of the convective core. Since this affects a very small 
part of the core of the model the effects on its global proper- 
ties are modest; however, it may have some influence on its 
pulsational properties, particularly for g modes or interface 
modes predominantly trapped near the boundary of the core, 
possibly offering the potential for asteroseismic diagnostics 
of this behaviour. Significant se nsitivity to the tre atment of 
this region was indeed found bv lMova et al.l (12007); a simi- 
lar model, but computed without such treatment of semicon- 

^ Owing to slight convergence problems in this calculation, V^d — 
Vgd is in fact positive at the point identified as the convective boundary. 
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vection (and with somewhat inadequate mesh resolution of 
the critical region) displayed far larger frequency differences 
than did the model considered in Fig. [T| between oscillation 
codes using different numerical techniques, particularly for 
modes with g-mode like behaviour and partly trapped near 
the edge of the core. 

The second aspect of semiconvection concerns the base 
of the convective envelope in models with diffu sion and set- 
tling o f helium and heavy elements. A s noted by lBahcall et al.^ 
( 2001 ) and discussed in more detail bv lChristen sen-Dalsgaard 
(120071) the gradient in the heavy-element abundance Z es- 
tablished by settling beneath the convection zone leads to 
convective instability, e.g., in IMq models of age somewhat 
higher than the present Sun. Complete mixing of the unsta- 
ble region removes the gradient and hence the cause of the 
instability, leading to an uncertain situation, characteristic of 
semiconvection. As above, a perhaps reasonable assumption 
might be that partial mixing takes place to establish com- 
position gradients that ensure neutral convective stability in 
the critical region, with Vjad = Vad. Since the opacity de- 
pends both on X and Z the mixing must consistently affect 
both abundances. I have attempted to implement this by in- 
cluding a turbulent diffusivity, obviously common to all ele- 
ments, determined iteratively as a function of depth beneath 
the convection zone such that the resulting profiles of X and 
Z lead to neutral stability. 



7 Concluding remarks 

No stellar evolution code is probably ever finished or fully 
tested. ASTEC has certainly proved useful in a number of 
applications, and the results for the Sun, as applied to helio- 
seismology, are perhaps reasonably reliable, at least within 
the framework of 'standard solar modelling'. It is obvious, 
however, that application to the increasingly accurate and 
detailed asteroseismic data will require further development. 
The tests provided through the ESTA collaboration and ex- 
tended through the HELAS Coordination Action are cer- 
tainly most valuable in this regard. 
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A Treatment of diffusion and settling 

For completeness I reproduce the expressions used to compute the 
diffusion and settling coef ficients. These are largely obtained from 
iMichaud and Proffit3 1 ll993f) . although with minor mo dific ations. 

In terms of the variable % introduced in Section [iTI Eqs [8l and [9l 
become 

dX, ■^M,^^^^_^^^^^ ^^^^ 



dx 9k 

For hydrogen = 1) 



logAHHe(0.7-h0.3X) \A 8 



9v^^(l-X-Z). 



(30) 



(31) 



where 

„ 15/2m„\'/^4^^ 

D = 



16 V 5;i: 



(32) 



niu being the atomic mass unit, Boltzmann's constant and e the elec- 
tron charge. Also, logAnHe is the Coulomb logarithm, approximated 
by 



log AHHe = - 19.95 - ^ log p - i log f j + ^ log r . 



(33) 



In the above expressions, cgs units are used and 'log' is natural loga- 
rithm. The diffusion coefficient is given by 



(3+X) 



plogAHHe(0.7-h0.3X) (\+X){-i + 5X) 
For a trace element k, with charge Zk and atomic mass Ak, 



1/2 



(34) 



(35) 



where 

'^k = XiA^mCkVi - AyHeCMe) . (36) 

Here Qh is t he collision integral between element k and hydrogen, 
which IMichau d and Proffitt ( 1993) approximate, fitting numerically 
determined values, as 

Cm = Y^log[exp(1.21ogA,H) + l] , (37) 
where 

AM=AHHe-logZ*+log2; (38) 

similarly the collision integral Qh between element k and helium is 
evaluated as 



Cme = Y^log[exp(1.21ogAtHe) + l]- 
where 

AA He=AHHe-l0gZj: . 

Also, 

AjAh 



A;,H = 



Ak+A-a 



(39) 



(40) 



(41) 



is the reduced mass in the element k, hydrogen system, Ah being the 
atomic mass of hydrogen; the reduced mass A^- He involving helium is 
defined similarly. 

As presented bv Mi chaud and Proffittl ( II993[) the diffusion velocity 
of trace elements depends on the gradient in the hydrogen abundance. 
To incorporate this in the formalism described by Eqs |29] and |30l I 
express dXjdr in terms of Y\ and write the coefficient ^ as 



Here 



/2^]2[l+Zk-Ak{5X + 3)/4] 



5i/2z2(<%+QHeA, 
0.54(4.75X-h2.25)^\ Gm 



1/2- 



(42) 



(43) 



log AHHe + 5 ) p 



2Zi, + 5(\+X) 



^1 5'/2pZf(^, +QHeA!(?J (1 +X){5X + 3) 



and 



{Anr^pf IBT^I^ 2Zk + 5{\+X) 



'^k 



oj, 5'/2pZf (l+X)(5X + 3) X(<r,+C,HeA:/?,; 



(44) 
-0.23 . 
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B Mixing-length formulation 



The calculation of V in convective regions is carried out us ing the 
[Vitense (1953) mixing length theory, in the form given by IGougl] 
n977i) . This is expressed in terms of 



iK/pc„)^Hp 



-1/2 



and 



1/2 



(45) 



(46) 



(47) 



Here Hp = p/{pg) is the pressure scale height, g = Gm/r^ is the grav- 
itational acceleration, 5 = — (51ogp/(9logr)p, K = /(3Kp) is 
the radiative conductivity and Cp is the specific heat at constant pres- 
sure. Also £ is the mixing length, which is taken to be a constant mul- 
tiple of the pressure scale height, £ = aui.Hp, and t] and <P are geo- 
metrical quantities, assumed be constant, which are related to the as- 
pect ratio of the conv ective cells; for t] = \/2/9 and <P = 2we. get the 
IBohm-Vitensd il958t) expressions (these values are used in the compu- 
tation). Then 



V-V,d=y(F-hA)(V„d-V,d) 
where Y is the positive root of 



1 



+ + AY -1 = . 



3AA 

The general solution to Ea.l49lis 



Y =A 



where 



--X{\-X)--X 
A ^ ' X 



x = A[Y+[f+X\\-Xf]'l^'^ 



1/3 



and 

r=x 



3 
2A2 



(48) 



(49) 



(50) 



(51) 



(52) 



For large and small A asymptotic expressions for Y are easily found. 
For A » 1 



\-A- 



and for A ^ 1 

y~(3AA)'/^^ 



2-\\a-' 
3A 



(53) 



(54) 



Eauation[53lis used when A > 15 and Ea.l54lwhen A < 10 ^. At these 
points the relative differences between tlie asymptotic values of Y and 
those found from Eq.[50]are less than 5 x 10^^. 



